

path <- "~/Dropbox/Friends/Yiqing and Haixiao/State_capacity/Replication"
setwd(path)
library(foreign)

#################################
## Figure 1
#################################

library(RColorBrewer) 
mycol <- c("gray50",1)
d<-read.dta("Data_full.dta");head(d)
par(mar=c(4.5,4.5,1,1))
plot(1,type="n",xlim=c(-4,4),ylim=c(0,0.5),xlab="",ylab="", cex.axis=1.5)
lines(density(d$Capacity[which(d$demo==0)]),col=mycol[1],lty=2,lwd=5)
lines(density(d$Capacity[which(d$demo==1)]),col=mycol[2],lwd=4)
mtext("State Capacity",1,line=3,cex=1.5)
mtext("Density",2,line=3,cex=1.5)
text(-2.2,0.4,"Autocracy",col=mycol[1],cex=1.5)
text(1.7,0.4,"Democracy",col=mycol[2],cex=1.5)
box()


## Raw dynmaics
d<-read.dta("fg_rawdyn.dta")
par(mar=c(4.5,4.5,1,1))
plot(1,type="n",axes=F,xlab="",ylab="",ylim=c(-0.45,0.2),xlim=c(-10,10))
abline(v=0,col="gray50",lty=2,lwd=4)
lines(d$demo_yr,d$Capacity,lwd=4, col = mycol[2])
points(d$demo_yr,d$Capacity,pch=16,cex=2, col = mycol[2])
axis(2,at=seq(-0.4,0.2,0.2), cex.axis=1.5) 
axis(1,at=seq(-8,8,4), cex.axis = 1.5)
mtext("Year(s) Since Democratization",1,line=3, cex = 1.5)
mtext("Average State Capacity",2,line=3, cex = 1.5)
box()

#################################
## Figure 2 Estimates
#################################

d <- cbind.data.frame("x0" = 1:7,
                      "coef" = c(24, 23, 23, 28, 26, 26, 22)/1000,
                      "se" = c(10,10,11,10,9,10,11)/1000)
d$ci1 <- d$coef - 1.96 * d$se
d$ci2 <- d$coef + 1.96 * d$se
labels <- c("None",
            "Population\nEthnic Frac.",
            "Income\nAgriculture\nLiteracy",
            "Minerals\nResource Rents",
            "Aid\nExports",
            "Interstate Rivalry\nTerritorial Threat\nWar",
            "All")
par(mar = c(1, 3.5, 0.2, 1))
plot(1, type = "n", xlim = c(0.5,7.5), ylim = c(-0.05, 0.08), axes = FALSE,
     xlab = "", ylab = "")
for (i in seq(0,8, 0.5)) {
    lines(c(i,i),c(-0.02,0.1), col = "gray90", lwd = 1.2)
}
##abline(v = seq(0,8,0.5), col = "gray90")
abline(h = seq(-0.02,0.08,0.01), col = "gray90", lwd = 1.2)
box()
mtext("Effect of Democracy on State Capacity", 2, line = 2.5, cex = 0.9)
axis(2, at = seq(-0.02,0.08, 0.02), cex.axis = 0.7)
abline(h = 0, col = "gray80", lwd = 4, lty = 2)
points(d$x0, d$coef, pch = 16, cex = 2, col = mycol[2])
for (i in 1:7) {lines(c(i,i), c(d$ci1[i],d$ci2[i]), lwd = 4, col = mycol[2])}
text(1.8,-0.02, "Control Variables", pos = 1, font = 2,
     col = "gray70", cex = 1.2)
axis(1, at= 1:7, labels = labels, cex.axis = 0.7, tick = FALSE, line = -2)

#############################################
## Figure 3. Contestation and Participation
#############################################

library(foreign)
library(RColorBrewer) 
mycol<-brewer.pal(3,"Set2")[c(2,1,3)]
d<-read.dta("Data_full.dta")
d<-d[,c("ccode","year","Capacity","contdim","partdim","demo")]
d<-na.omit(d)
dem<-d[which(d$demo==1),]
aut<-d[which(d$demo==0),]
pcex<-0.6; tcex<-1.5;

## plot
X<-"contdim"
par(lend=1,mar=c(4.5,4.5,1,1))
plot(1,type="n",xlim=c(0,1),ylim=c(-3,3),xlab="",ylab="", cex.axis = 1.5)
points(aut[,X],aut$Capacity,col=mycol[1],pch=1,cex=pcex)
points(dem[,X],dem$Capacity,col=mycol[2],pch=2,cex=pcex)
x0<-seq(0,0.99,by=0.01)
fit <- loess(d$Capacity ~ d[,X], span=0.6)
prd <- predict(fit, x0, se=TRUE)
x.poly <- c(rev(x0), x0)
y.poly <- c(rev(prd$fit+1.96*prd$se.fit), (prd$fit-1.96*prd$se.fit))
polygon(x.poly, y.poly,col="#66666670", border=NA)
lines(x0,prd$fit,lwd=3)
mtext("State Capacity",2,3,cex=tcex)
mtext("Contestation",1,3,cex=tcex)
legend("topleft",c("Autocracy","Democracy"),cex=1.5,pch=c(1,2),
       col=mycol,bty="n")

X<-"partdim"
par(lend=1,mar=c(4.5,4.5,1,1))
plot(1,type="n",xlim=c(0,1),ylim=c(-3,3),xlab="",ylab="", cex.axis = 1.5)
points(aut[,X],aut$Capacity,col=mycol[1],pch=1,cex=pcex)
points(dem[,X],dem$Capacity,col=mycol[2],pch=2,cex=pcex)
x0<-seq(0,0.99,by=0.01)
fit <- loess(d$Capacity ~ d[,X], span=0.6)
prd <- predict(fit, x0, se=TRUE)
x.poly <- c(rev(x0), x0)
y.poly <- c(rev(prd$fit+1.96*prd$se.fit), (prd$fit-1.96*prd$se.fit))
polygon(x.poly, y.poly,col="#66666670", border=NA)
lines(x0,prd$fit,lwd=4)
mtext("State Capacity",2,3,cex=tcex)
mtext("Participation",1,3,cex=tcex)
legend("topleft",c("Autocracy","Democracy"),cex=1.5,pch=c(1,2),
       col=mycol,bty="n")


######################## Appendix #########################

##################################
## Figure A1. Scatterplot Matrix
##################################

library(foreign)
library(car)
d<-read.dta("fg_scatterplot.dta");names(d)

scatterplotMatrix(~Capacity+rpc+pubcorr+logmilper,pch=16,cex=0.6,col=c(1,2,"#99999950"),lwd=3,data=d,main="")

scatterplotMatrix(~Capacity_res+rpc_res+pubcorr_res+logmilper_res,pch=16,cex=0.6,col=c(1,2,"#99999950"),lwd=3,data=d,main="")

########################################################################
## Figure A2.  Democratization (number of democracies and autocracies)
#########################################################################

d<-read.dta("fg_demo.dta");head(d)
par(mar=c(3.5,3.5,1,1))
plot(d$year[d$demo==0],d$countrynum[d$demo==0],type="l",ylim=c(0,120),xlab="",ylab="",col="gray70",lty=2,lwd=5)
lines(d$year[d$demo==1],d$countrynum[d$demo==1],lwd=3)
mtext("Year",1,line=2,cex=1.2)
mtext("Number of countries",2,line=2,cex=1.2)
text(1990,100,"Autocracies",col="gray50",cex=1.5)
text(1985,35,"Democracies",cex=1.5)

########################################################################
## Figure A3. Year of democratization / reversal
########################################################################

d<-read.dta("fg_demoyr.dta");head(d)
# total counts
t<-data.frame(table(d$year));colnames(t)<-c("year","counts")
# democratization counts
s<-data.frame(table(d$year[d$check==1]));colnames(s)<-c("year","counts_demo")
# merge to a data frame without missing years
year<-c(1960:2010)
m<-data.frame(year)
m<-merge(m,t,by="year",all=T)
m<-merge(m,s,by="year",all=T)
m$counts[is.na(m$counts)==1]<-0
m$counts_demo[is.na(m$counts_demo)==1]<-0
# plotting
par(mar=c(3.5,3.5,1,1))
barplot(m$counts,space=c(0.7),border=1,ylim=c(0,8),col="white",names.arg=m$year)
barplot(m$counts_demo,space=c(0.7),border=1,add=TRUE,col="gray40")
axis(1,at=seq(1,90,8.5),labels=F)
box();mtext("Number of Regime Changes",2,2,cex=1.2);mtext("Year",1,2,cex=1.2)
legend("topright",c("Democratization","Reversal"),cex=1.2, fill=c("gray40","white"),seg.len=1,bty="n")


#########################
## Figure A4. maps
#########################

##install.packages("rworldmap")
library(rworldmap)
library(foreign)
d<-read.dta("fg_map.dta")
head(d)
d<-rbind(d,c("S. Sudan","SSD",1970,0),c("S. Sudan","SSD",1985,0),c("S. Sudan","SSD",2000,0)) 
d<-rbind(d,c("Russia","RUS",1970,0),c("Russia","RUS",1985,0))
d<-rbind(d,c("Ukraine","UKR",1970,0),c("Ukraine","UKR",1985,0))
d<-rbind(d,c("Kazakhstan","KAZ",1970,0),c("Kazakhstan","KAZ",1985,0))
d<-rbind(d,c("Belars","BUR",1970,0),c("Belars","BUR",1985,0))
d<-rbind(d,c("Germany","DEU",1970,1),c("Germany","DEU",1985,1),c("Germany","DEU",2000,1))
d<-rbind(d,c("","BLR",1970,0),c("","BLR",1985,0))
d<-rbind(d,c("","TKM",1970,0),c("","TKM",1985,0))
d<-rbind(d,c("","AZE",1970,0),c("","AZE",1985,0))
d<-rbind(d,c("","UZB",1970,0),c("","UZB",1985,0))
d<-rbind(d,c("","KGZ",1970,0),c("","KGZ",1985,0))
d<-rbind(d,c("","TJK",1970,0),c("","TJK",1985,0))
d<-rbind(d,c("","NAM",1970,0),c("","NAM",1985,0))
d<-rbind(d,c("","LVA",1970,0),c("","LVA",1985,0))
d<-rbind(d,c("","LTU",1970,0),c("","LTU",1985,0))
d<-rbind(d,c("","ROU",1970,0),c("","ROU",1985,0))
d<-rbind(d,c("","CZE",1970,0),c("","CZE",1985,0))
d<-rbind(d,c("","SVK",1970,0),c("","SVK",1985,0))
d<-rbind(d,c("","MDA",1970,0),c("","MDA",1985,0))
d<-rbind(d,c("","BGD",1970,0))
d<-rbind(d,c("","AGO",1970,0))
d<-rbind(d,c("","MOZ",1970,0))
s70<-d[which(d$year==1970),]  # 131 countries
s85<-d[which(d$year==1985),]  # 141
s00<-d[which(d$year==2000),]  # 160
map.demo.70<-joinCountryData2Map(s70,joinCode="ISO3",nameJoinColum="isocode")
map.demo.85<-joinCountryData2Map(s85,joinCode="ISO3",nameJoinColum="isocode")
map.demo.00<-joinCountryData2Map(s00,joinCode="ISO3",nameJoinColum="isocode")
##color
mycol<-c("gray30","gray80")

## 1970
par(mar=c(0,1,1,1))
mapCountryData(map.demo.70,nameColumnToPlot="demo",catMetho="categorical",mapTitle="",aspect=1,addLegend=F,
               colourPalette=mycol,borderCol="gray10",ylim=c(-20,55))
legend("top",c("Democracy","Autocracy"),fill=mycol[c(2,1)],bty="n",ncol=2);box()

## 1985
par(mar=c(0,1,1,1))
mapCountryData(map.demo.85,nameColumnToPlot="demo",catMetho="categorical",mapTitle="",aspect=1,addLegend=F,
               colourPalette=mycol,borderCol="gray10",ylim=c(-20,55))
legend("top",c("Democracy","Autocracy"),fill=mycol[c(2,1)],bty="n",ncol=2);box()

## 2000
par(mar=c(0,1,1,1))
mapCountryData(map.demo.00,nameColumnToPlot="demo",catMetho="categorical",mapTitle="",aspect=1,addLegend=F,
               colourPalette=mycol,borderCol="gray10",ylim=c(-20,55))
legend("top",c("Democracy","Autocracy"),fill=mycol[c(2,1)],bty="n",ncol=2);box()


########################################################################
## Figure A5. dynamic effect
########################################################################

# dynamic effect
(d<-read.csv("fg_dyn.csv"))
d<-d[which(d$period>=-9&d$period<=10),] # 10 periods before; 10 periods after
d$CI1_fe<-d$coef_fe-1.96*d$se_fe
d$CI2_fe<-d$coef_fe+1.96*d$se_fe
d$CI1_lag2<-d$coef_lag2-1.96*d$se_lag2
d$CI2_lag2<-d$coef_lag2+1.96*d$se_lag2
x0<-d$period

# fixed effects
par(mar=c(3,3,1,1))
plot(1,type="n",xlim=c(-9.5,10.5),ylim=c(-0.3,0.42),xlab="",ylab="",axes=F);box()
abline(h=0,col="gray",lwd=2)
abline(v=0,col=2,lty=2)
polygon(c(rev(x0),x0),c(rev(d$CI2_fe),d$CI1_fe),col="#99999980",border=NA)
lines(x0,d$coef_fe,lwd=5)
mtext("Year(s) since Democratization",1,line=2,cex=1.2);axis(1,at=c(seq(-10,10,by=2),1));axis(2)
mtext("Coefficeint",2,line=2,cex=1.2)

# plus lages
par(mar=c(3,3,1,1))
plot(1,type="n",xlim=c(-9.5,10.5),ylim=c(-0.3,0.42),xlab="",ylab="",axes=F);box()
abline(h=0,col="gray",lwd=2)
abline(v=0,col=2,lty=2)
polygon(c(rev(x0),x0),c(rev(d$CI2_lag2),d$CI1_lag2),col="#99999980",border=NA)
lines(d$period,d$coef_lag2,lwd=5)
mtext("Year(s) since Democratization",1,line=2,cex=1.2);axis(1,at=c(seq(-10,10,by=2),1));axis(2)
mtext("Coefficeint",2,line=2,cex=1.2)


############################
## Figure A6. Cases (Peru)
############################


d<-read.dta("fg_modelfit.dta")
head(d)
s<-d[which(d$ccode==135),]
years<-unique(s$year)
par(mar=c(2,3,1,1),mfcol=c(2,1))
plot(1,type="n",xlab="",ylab="",xlim=c(1960,2010),ylim=c(-1.05,0.95),axes=F);box();axis(1)
axis(2,at=seq(-0.8,0.8,0.4));mtext("Capacity",2,2)
for (i in years) {if (s$demo[which(s$year==i)]==1) {polygon(c(i+1,i,i,i+1),c(2,2,-2,-2),col="#AAAAAA50",border=NA)}}
lines(s$year,s$Cap_fe,col="gray50",lwd=3,lty=2)
lines(s$year,s$Cap_td1,col="gray50",lwd=3,lty=3)
lines(s$year,s$Cap_td2,col="gray50",lwd=3,lty=4)
lines(s$year,s$Cap_td3,col="gray50",lwd=3,lty=5)
lines(s$year,s$Capacity,lwd=3)
text(1959,-0.65,"FE",pos=4,col="gray40")
text(1965,-0.7,"FE+linear trend",pos=4,col="gray40")
text(1993,0.6,"FE+quadratic/cubic trend",pos=4,col="gray40")
text(1982,0.55,"Actual",pos=1,col=1)
# lags
plot(1,type="n",xlab="",ylab="",xlim=c(1960,2010),ylim=c(-1.05,0.95),axes=F);box();axis(1)
axis(2,at=seq(-0.8,0.8,0.4));mtext("Capacity",2,2)
for (i in years) {
    if (s$demo[which(s$year==i)]==1) {
        polygon(c(i+1,i,i,i+1),c(2,2,-2,-2),col="#AAAAAA50",border=NA)
    }
}
lines(s$year,s$Cap_lag2,lwd=3,col=2,lty=2)
lines(s$year,s$Capacity,lwd=3)
text(1988.5,0.3,"FE+2 lags",pos=1,col=2)
text(1982,0.55,"Actual",pos=1,col=1)

#############################
# Figure A7. IV first stage 
#############################

library(foreign)
d<-read.dta("fg_firststage.dta")
par(lend=1,mar=c(3.5,3.5,1,1))
plot(1,type="n",xlim=c(0,1),ylim=c(0,1),xlab="",ylab="")
rug(d$demo_wave[which(d$demo==0)],col=1,side=1)
rug(d$demo_wave[which(d$demo==1)],col=1,side=3)
# lowess with confidence interval
x0<-seq(0,1,by=0.02)
fit <- loess(d$demo ~ d$demo_wave, span=0.6)
prd <- predict(fit, x0, se=TRUE)
x.poly <- c(rev(x0), x0)
y.poly <- c(rev(prd$fit+1.96*prd$se.fit), (prd$fit-1.96*prd$se.fit))
polygon(x.poly, y.poly,col="#55555550", border=NA)
lines(x0,prd$fit,lwd=4)
mtext("Democracy",2,2,cex=1.5)
mtext("Democractic Wave",1,2,cex=1.5)

par(lend=1,mar=c(3.5,3.5,1,1))
plot(1,type="n",xlim=c(-0.25,0.25),ylim=c(-0.45,0.45),xlab="",ylab="")
points(d$demo_wave_res,d$demo_res,col="#99999950",pch=1,cex=1)
# lowess with confidence interval
x0<-seq(-0.25,0.25,by=0.01)
fit <- loess(d$demo_res ~ d$demo_wave_res, span=0.6)
prd <- predict(fit, x0, se=TRUE)
x.poly <- c(rev(x0), x0)
y.poly <- c(rev(prd$fit+1.96*prd$se.fit), (prd$fit-1.96*prd$se.fit))
polygon(x.poly, y.poly,col="#55555550", border=NA)
lines(x0,prd$fit,lwd=4)
mtext("Residualized Democracy",2,2,cex=1.5)
mtext("Residualized Democractic Wave ",1,2,cex=1.5)


##################################
## Figure A8. Mediation Analysis
##################################

## ##install.packages("mediation")
## library(mediation)
## library(foreign)

## ## prepare
## d<-read.dta("Data_mediation.dta")
## names(d)

## d$ccode<-as.factor(d$ccode)
## d$year<-as.factor(d$year)
## cc<-model.matrix(~ccode-1,data=d)[,-1]
## yy<-model.matrix(~year-1,data=d)[,-1]
## d<-cbind(d,cc,yy)

## ## variable names
## Y<-"Capacity"
## D<-"demo_l1"
## M1<-"contdim_l1"
## M2<-"partdim_l1"
## cluster<-d[,"ccode"]
## (covars<-c(colnames(cc),colnames(yy),"Capacity_l1","Capacity_l2",names(d)[7:18]))
## sims<-1000

## ## contestation
## mod.m1<-lm(paste(M1,"~",do.call(paste,c(as.list(c(D,covars)),sep=" + "))),data=d)
## mod.y1<-lm(paste(Y,"~",do.call(paste,c(as.list(c(D,M1,covars)),sep=" + "))),data=d)
## med1<-mediate(mod.m1,mod.y1,sims=sims,treat=D,mediator=M1,robustSE=TRUE,
##               conf.level = 0.90)
## summary(med1)
## plot(med1)

## ## participation
## mod.m2<-lm(paste(M2,"~",do.call(paste,c(as.list(c(D,covars)),sep=" + "))),data=d)
## mod.y2<-lm(paste(Y,"~",do.call(paste,c(as.list(c(D,M2,covars)),sep=" + "))),data=d)
## med2<-mediate(mod.m2,mod.y2,sims=sims,treat=D,mediator=M2,robustSE=TRUE,
##               conf.level = 0.90)
## summary(med2)
## plot(med2)

## ## sensitivity analysis
## sens.out1 <- medsens(med1, rho.by = 0.1, sims = 500, effect.type = c("both"))
## sens.out2 <- medsens(med2, rho.by = 0.1, sims = 500, effect.type = c("both"))

## ## save results
## save(d, med1, med2, sens.out1, sens.out2, file="result_mediation.RData")


## plotting
library(mediation)
load("fg_mediation.RData")

pcex=1.8;s=0.05
par(mar=c(4.5,6.5,1,0.5),las=1)
plot(1,xlim=c(-0.04,0.06),ylim=c(0.5,3.6),type="n",
     axes=F,xlab="",ylab="");box();axis(1, cex.axis = 1.5);
abline(v=0,col="gray50",lty=2,lwd=2)
lines(med1$tau.ci,c(1,1),lwd=3);points(med1$tau.coef,1,pch=16,cex=pcex)
lines(med1$z0.ci,c(2,2)-s,lwd=3);points(med1$z0,2-s,pch=16,cex=pcex)
lines(med1$d0.ci,c(3,3)-s,lwd=3);points(med1$d0,3-s,pch=16,cex=pcex)
lines(med2$z0.ci,c(2,2)+s,lwd=3,col="gray50",lty=2);points(med2$z0,2+s,pch=17,cex=pcex,col="gray50")
lines(med2$d0.ci,c(3,3)+s,lwd=3,col="gray50",lty=2);points(med2$d0,3+s,pch=17,cex=pcex,col="gray50")
mtext("The Effect of Democarcy on State Capacity",1,3,cex=1.5)
axis(2,at=1:3,labels=c("Total\nEffect","Direct\nEffect","Mediation\nEffect"), cex.axis = 1.5)
legend("topright",c("Contestation","Participation"),cex=1.5,pch=c(16,17),
       col=c(1,"gray50"),lwd=3,seg.len=3,lty=c(1,2),bty="n")


par(mfcol=c(1,2))
plot(sens.out1, ask = FALSE, xlim= c(-0.1,0.1), ylim = c(-0.15, 0.15))

par(mfcol=c(1,2))
plot(sens.out2, ask = FALSE, xlim= c(-0.1,0.1), ylim = c(-0.15, 0.15))




